#loading packages
library('FactoMineR')
library('factoextra')
library("gplots")
library("ggrepel")
library("dplyr")
library("tidyverse")
library("ggplot2")
library("ggthemes")
library("forcats")
library("reshape2")

#datasets

#Chapter 2
load(file="Mapping Religious Accommodation in the EU/ReligiousAccommodationEurope.rdata") #Dataset used for the MCA - Country MCA

load(file="Mapping Religious Accommodation in the EU/accommeans.rdata") #Dataset used for the Score analysis - Country MCA

#MCA for Countries - ReligiousAccomodationEurope dataset

##load dataset
load(file="Mapping Religious Accommodation in the EU/ReligiousAccommodationEurope.rdata") #Dataset used for the MCA - Country MCA

##MCA
mca <- MCA(Accom_active, graph = FALSE)
print(mca)

##Variable Coordinates Contribution
var <- get_mca_var(mca)
dt <-var$coord
dt <- data.frame(dt[1:36, 1:2])

##Adding Level of Restriction and Contribution
Restriction <- c(0, 1, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 0, 2, 3, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 0, 2, 0, 2, 0, 1, 2, 3)
dt$Restrictions <- Restriction
dt$Mean <- round(rowMeans(dt[,c('Dim.1', 'Dim.2')], na.rm=TRUE), 2)
dt$Size <- abs(dt$Mean)
size <- c(0, 0, 0, 3, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 2, 0, 0, 0, 1, 0, 0, 1, 1)
dt$size <- size

rownames(dt) <- c("Rites of Passage", "Rites of Passage ", "Diet Laws", "Diet Laws ", "Diet Laws  ", "Burial", "Burial ", "Burial  ", "Rel. Clothing", "Rel. Clothing ", "Rel. Clothing  ", "Rel. Symbols", "Rel. Symbols ", "Public Holidays", "Public Holidays ", "Public Holidays  ", "Private Holidays", "Private Holidays ", "Clergy Jails", "Clergy Jails ", "Clergy Jails  ", "Clergy Military Base", "Clergy Military Base ", "Clergy Military Base  ", "Clergy Military Base   ", "Clergy Hospitals", "Clergy Hospitals ", "Clergy Hospitals  ", "Public Preaching", "Public Preaching ", "Proselytizing", "Proselytizing ", "Religious Ed. Mandatory", "Religious Ed. Mandatory ", "Religious Ed. Mandatory  ", "Religious Ed. Mandatory   ")

##Visualising Variable Categories

ggplot(dt, aes(x=Dim.1, y=Dim.2)) +
  geom_point(aes(color = factor(Restrictions), size=Size)) +
  labs(x = "Privatisation of Religion", y = "Regulation of Religion") +
  geom_hline(yintercept=0, color="grey") +
  geom_vline(xintercept=0, color="grey") +
  geom_text_repel(label=rownames(dt), size=5, family = "Times New Roman")+
  theme_classic(base_family = "Times New Roman") +
  labs(title = "Figure 1 Religious Accommodation Variable Categories", subtitle="Contribution of each variable where the dot size represents their importance for each dimension.", caption= "Source: Religion and the State Round 3, Global Restrictions on Religion. Own Elaboration") +
  theme(title= element_text(size = 15), text= element_text(size= 8, family = "Times New Roman"), legend.position = c(0.2, 0.25), legend.justification = c("right", "top"),  legend.text = element_text(size=15), axis.text=element_text(size=12), axis.title=element_text(size=14)) +
  guides(color = guide_legend(override.aes = list(size=5))) +
  labs(color="Degree of Intervention") +
  scale_color_manual(labels = c("No", "Little", "Somewhat", "High"), values=c("#488f31", "#9fc08f", "#fda0c3", "#fa3497"))+
  guides(size=FALSE)

##Extracting Individuals Data for Country Categorisation
df <-mca$ind$coord
df <- data.frame(df[1:216, 1:2])

##Selecting and summarising per countries
df_test<- data.frame(df[1:8, 1:2])
df_test <- summarise_each(df_test, funs(mean))
rownames(df_test) <- c("Austria")

Belgium <- data.frame(df[9:16, 1:2])
Belgium <- summarise_each(Belgium, funs(mean))
rownames(Belgium) <- c("Belgium")
Clusters <- dplyr::union(df_test, Belgium)

Bulgaria <- data.frame(df[17:24, 1:2])
Bulgaria <-summarise_each(Bulgaria, funs(mean))
rownames(Bulgaria) <- c("Bulgaria")
Clusters <- dplyr::union(Clusters, Bulgaria)

Croatia <- data.frame(df[25:32, 1:2])
Croatia <- summarise_each(Croatia, funs(mean))
rownames(Croatia) <- c("Croatia")
Clusters <- dplyr::union(Clusters, Croatia)

Czech <- data.frame(df[33:40, 1:2])
Czech <- summarise_each(Czech, funs(mean))
rownames(Czech) <- c("Czech")
Clusters <- dplyr::union(Clusters, Czech)

Denmark <- data.frame(df[41:48, 1:2])
Denmark <- summarise_each(Denmark, funs(mean))
rownames(Denmark) <- c("Denmark")
Clusters <- dplyr::union(Clusters, Denmark)

Estonia <- data.frame(df[49:56, 1:2])
Estonia <- summarise_each(Estonia, funs(mean))
rownames(Estonia) <- c("Estonia")
Clusters <- dplyr::union(Clusters, Estonia)

Finland <- data.frame(df[57:64, 1:2])
Finland <- summarise_each(Finland, funs(mean))
rownames(Finland) <- c("Finland")
Clusters <- dplyr::union(Clusters, Finland)

France <- data.frame(df[65:72, 1:2])
France <- summarise_each(France, funs(mean))
rownames(France) <- c("France")
Clusters <- dplyr::union(Clusters, France)

Germany <- data.frame(df[73:80, 1:2])
Germany <- summarise_each(Germany, funs(mean))
rownames(Germany) <- c("Germany")
Clusters <- dplyr::union(Clusters, Germany)

Greece <- data.frame(df[81:88, 1:2])
Greece <- summarise_each(Greece, funs(mean))
rownames(Greece) <- c("Greece")
Clusters <- dplyr::union(Clusters, Greece)

Hungary <- data.frame(df[89:96, 1:2])
Hungary <- summarise_each(Hungary, funs(mean))
rownames(Hungary) <- c("Hungary")
Clusters <- dplyr::union(Clusters, Hungary)

Ireland <- data.frame(df[97:104, 1:2])
Ireland <- summarise_each(Ireland, funs(mean))
rownames(Ireland) <- c("Ireland")
Clusters <- dplyr::union(Clusters, Ireland)

Italy <- data.frame(df[105:112, 1:2])
Italy <- summarise_each(Italy, funs(mean))
rownames(Italy) <- c("Italy")
Clusters <- dplyr::union(Clusters, Italy)

Latvia <- data.frame(df[113:120, 1:2])
Latvia <- summarise_each(Latvia, funs(mean))
rownames(Latvia) <- c("Latvia")
Clusters <- dplyr::union(Clusters, Latvia)

Lithuania <- data.frame(df[121:128, 1:2])
Lithuania <- summarise_each(Lithuania, funs(mean))
rownames(Lithuania) <- c("Lithuania")
Clusters <- dplyr::union(Clusters, Lithuania)

Luxembourg <- data.frame(df[129:136, 1:2])
Luxembourg <- summarise_each(Luxembourg, funs(mean))
rownames(Luxembourg) <- c("Luxembourg")
Clusters <- dplyr::union(Clusters, Luxembourg)

Malta <- data.frame(df[137:144, 1:2])
Malta <- summarise_each(Malta, funs(mean))
rownames(Malta) <- c("Malta")
Clusters <- dplyr::union(Clusters, Malta)

Netherlands <- data.frame(df[145:152, 1:2])
Netherlands <- summarise_each(Netherlands, funs(mean))
rownames(Netherlands) <- c("Netherlands")
Clusters <- dplyr::union(Clusters, Netherlands)

Poland <- data.frame(df[153:160, 1:2])
Poland <- summarise_each(Poland, funs(mean))
rownames(Poland) <- c("Poland")
Clusters <- dplyr::union(Clusters, Poland)

Portugal <- data.frame(df[161:168, 1:2])
Portugal <- summarise_each(Portugal, funs(mean))
rownames(Portugal) <- c("Portugal")
Clusters <- dplyr::union(Clusters, Portugal)

Romania <- data.frame(df[169:176, 1:2])
Romania <- summarise_each(Romania, funs(mean))
rownames(Romania) <- c("Romania")
Clusters <- dplyr::union(Clusters, Romania)

Slovakia <- data.frame(df[177:184, 1:2])
Slovakia <- summarise_each(Slovakia, funs(mean))
rownames(Slovakia) <- c("Slovakia")
Clusters <- dplyr::union(Clusters, Slovakia)

Slovenia <- data.frame(df[185:192, 1:2])
Slovenia <- summarise_each(Slovenia, funs(mean))
rownames(Slovenia) <- c("Slovenia")
Clusters <- dplyr::union(Clusters, Slovenia)

Spain <- data.frame(df[193:200, 1:2])
Spain <- summarise_each(Spain, funs(mean))
rownames(Spain) <- c("Spain")
Clusters <- dplyr::union(Clusters, Spain)

Sweden <- data.frame(df[201:208, 1:2])
Sweden <- summarise_each(Sweden, funs(mean))
rownames(Sweden) <- c("Sweden")
Clusters <- dplyr::union(Clusters, Sweden)

UK <- data.frame(df[209:216, 1:2])
UK <- summarise_each(UK, funs(mean))
rownames(UK) <- c("UK")
Clusters <- dplyr::union(Clusters, UK)

##Visualising Country Placements Biplot

ggplot(Clusters, aes(x=Dim.1, y=Dim.2)) +
  geom_point(color="#e43085", size=3) +
  labs(x = "Privatisation of Religion", y = "Regulation of Religion") +
  geom_hline(yintercept=0, color="grey") +
  geom_vline(xintercept=0, color="grey") + 
  geom_text_repel(label=rownames(Clusters), size=5, family = "Times New Roman") +
  theme_classic(base_family = "Times New Roman") +
  labs(title = "Figure 3 Religious Accommodation in the European Union", subtitle = "Average country placement according to Privatisation of Religion and Regulation of Religion.", caption= "Source: Religion and the State Round 3, Global Restrictions on Religion. Own Elaboration") +
  theme(title= element_text(size = 15), text= element_text(size= 8), axis.text = element_text(size=12), axis.title = element_text(size=14))

##Adding all values for ranking
load("Mapping Religious Accommodation in the EU/accommeans.rdata")
rel <- round(rel, digits=2)

names(rel) <- c("Rites of Passage","Dietary Laws", "Burial","Religious Clothing", "Religious Symbols","Public Holidays", "Private Holidays", "Clergy Jails","Clergy Military Base", "Clergy Hospitals", "Public Preaching","Proselytizing","Mandatory Religious Education")

rel <- rel[,c("Rites of Passage","Dietary Laws", "Burial","Religious Clothing", "Religious Symbols","Public Holidays", "Private Holidays", "Clergy Jails","Clergy Military Base", "Clergy Hospitals", "Public Preaching","Proselytizing","Mandatory Religious Education")]

rel$Sum <- rowSums(rel)
rel <- arrange(rel, desc(Sum))
rel = subset(rel, select = -c(Sum) )

rel$Country <- rownames(rel)

rel <- rel[,c("Country",names(sort(colSums(rel[,c(1:ncol(rel)-1)]), decreasing = T)))]

ord <- rel$Country
ord <- ord[length(ord):1]

rel1 <- melt(rel)

rel1 <- rel1 %>%
  mutate(Country = fct_relevel(Country, ord))

rel1$value <- round(rel1$value,1)

##Heatmap
rel1 %>%
  ggplot( aes(variable, Country, fill= value))+
  geom_tile(color = "white",
            lwd = 0.5,
            linetype = 1) +
  geom_text(aes(label = value), color = "white", size = 5, family = "Times New Roman")+
  xlab("Variable") +
  ylab("Country") +
  ggtitle("Figure 4 Average Country Scores in Original Datasets")+
  labs(fill = "Restrictions", subtitle="Brighter shades represent high levels of restrictions, while softer shades represent little to no restrictions") +
  scale_fill_gradient(low="white",high="#fa3497")+
  theme_classic(base_family = "Times New Roman") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 60, vjust = 0, hjust=0, size=12),
        plot.margin=unit(c(0.5,1,0.5,0),"cm"),
        plot.background = element_rect(color=NA), title= element_text(size = 13), text= element_text(size= 13), axis.text.y = element_text(size=12))+
  xlab("")+
  ylab("")+
  scale_x_discrete(position = "top")+
  labs(caption= "Source: Religion and the State Round 3, Global Restrictions on Religion. Own Elaboration")